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ABSTRACT 



Context. Supersonic turbulence is ubiquitous in the interstellar medium and plays an important role in contemporary star formation. 
Aims. To perform a high-resolution numerical simulation of supersonic isothermal turbulence driven by compressive large-scale 
forcing and to analyse various statistical properties. 

Methods. The compressible Euler equations with an external stochastic force iield dominated by rotation-free modes are solved with 
the piecewise parabolic method. Both a static grid and adaptive mesh refinement is used with an effective resolution A' = 768'. 
Results. After a transient phase dominated by shocks, turbulence evolves into a steady state with root mean square Mach number 
X 2.2 . . . 2.5, in which cloud-like structures of over-dense gas are surrounded by highly rarefied gas. The index of the turbulence energy 
spectrum function j8 « 2.0 in the shock-dominated phase. As the flow approaches statistical equilibrium, the spectrum flattens, with 
yS K 1.9. For the scaling exponent of the root mean square velocity fluctuation, we obtain y « 0.43 from the velocity structure functions 
of second order These results are well within the range of observed scaling properties for the velocity dispersion in molecular clouds. 
Calculating structure functions of order p = 1, ... ,5, we find for all scaling exponents significant deviations from the Kolmogorov- 
Burgers model proposed by Boldyrev. Our results are very well described by a general log-Poisson model with a higher degree of 
intermittency, which implies an influence of the forcing on the scaling properties. The spectral index of the quadratic logarithmic 
density fluctuation ps « 1.8. Contrary to previous numerical results for isothermal turbulence, we obtain a skewed probability density 
function of the mass density fluctuations that is not consistent with log-normal statistics and entails a substantially higher fraction 
of mass in the density peaks than implied by the Padoan-Nordlund relation between the variance of the density fluctuations and the 
Mach number 

Conclusions. Even putting aside further complexity due to magnetic fields, gravity or thermal processes, we question the notion that 
Larson-type relations are a consequence of universal supersonic turbulence scaling. For a genuine understanding, it seems necessary 
to account for the production mechanism of turbulence in the ISM. 

Key words. Hydrodynamics - Turbulence - Methods: numerical - ISM: kinematics and dynamics 



1. Introduction 



The question of the origin of supersonic inte r stellar 
turbulen ce remains largely op en (Elmegreen & Scalo 2004|; 
iBallester os-Paredes et al.l l2007h . For the numerical simulation 



Today it is agreed that supersonic turbulence occurs over a 
large range of length scales in the interstellar medium (ISM) 
and, p articularly, within molecular clouds (Elmegreen & Scalo,; 
12004 . This is inferred from different observations of molec- 
ular line broadening which indicate, firstly, a velocity disper- 
sion comparable to or greater than the speed of sound and, 
secondly, a power-law dependence of the velocity dispersion, 
Av oc i"^, in the range of length scales 0.1 p c < £ < 10 pc 
of the associated struc t ures (Larson|[T9 81: Falg arone et al.ll 19921 
iMiesch & Ballvl ll994': 'Bru nt & Heyeill200 2)."This has been in- 
terpreted as a man ifestatio n of t he universality of interstellar 
turbulence dHever & Brun^ [2004 ^ and bears important conse- 
quences on the theory of turbulence-regulated star formation, 
where Larson-type relations for the velocity dispersion are ex- 
plained as a direct consequence of the self-similarity of the 
turbulence cascade, with possible modifications of the scaling 
expone nt due to magnetic fields, self-gravity or thermal pro- 
cesses jP adoan & Nordlundl 12002 !: 'Mac Low & KlessenI l2004t 
iBallestero s-Pare des et alJl2007l; lMcKee & Ostrikeiil2007h . 



of interstellar turbulence, there is also the difficulty that en- 
ergy injection into the interstellar medium and processes lead- 
ing to star formation may occur on vastly disparate scales, 
which cannot be encompassed by a single numerical sim- 
ulation. This is why two major approaches to the problem 
have emerged: Either the generation of turbulence by in- 
stabilities at relatively large scales is computed or statisti- 
cally isotropic turbulence is simulated in a periodic box as 
a simple model for the small-scale regime. Typical examples 
for the first approach ar e the simulations of col li ding flows 
by Heitschetal. (200 6i). IVazquez-Semadeni et alj (l2007h and 
iHennebelle & Audill (I2007ir Recent simulations of driven turbu- 
len ce in a periodic box w ere presented by Pado an et al.l ( |2007l) 
and lKritsuk et al.l (l2007al) . While colliding flow simulations fol- 
low the idea of controlling the flow evolution in the compu- 
tational domain via initial and boundary conditions only, the 
meaning of turbulence simulations with a driving force are im- 
plicitly based on the notion of an upper numerical cutoflF scale. 
This means that one imagines turbulence being produced on 
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scales larger than the size of the computational box and the driv- 
ing force is intended to mimic the interactions with the larger 
scales in a heuristic fashion. Indeed, decomposing the hydrody- 
namic equations into a large-scale part and a small-scale part 
entails coupling via an energy flux that, in the most simplified 
manner, can be modelled as a random volume force acting upon 
the small-scale part Frisch ( 1995). Such a volume force can pro- 
duce shear or compression by exciting solenoidal (divergence- 
free) or dilatational (rotation-free) modes of the velocity field. 

In this article, we follow a hybrid approach. A statistically 
isotropic stochastic force field that excites mostly dilatational 
modes is applied in our simulations. Thereby, converging su- 
personic flows are produced randomly at large scales. From the 
collisions of these flows in combination with a small solenoidal 
component of the stochastic force, turbulence is generated and 
the flow evolves towards a steady state, in which energy injection 
due to the forcing is balanced by numerical dissipation. We refer 
to this scenario as compressively driven turbulence. Apart from 
the dominance of the rotation-free component, our method dif- 
fers in another aspec t form the forcing applied, for instance, by 
iKritsuk et al.1 (l2007al) and in many earlier simulations. Whereas 
they imprint a certain large-scale structure by a steady random 
force that is derived from some initial velocity field and merely 
replenishes the energy of the flow, we start with isothermal gas 
of uniform density at rest and let turbulence develop dynami- 
cally. The spatiotemporal correlation properties of the stochastic 
force field guarantee that large-scale structures evolve on a time 
scale that is in accordance with their characteristic velocity and 
size. Moreover, a statistically isotropic state is approached re- 
gardless of the initial condition, because temporal correlations 
decay exponentially. There is a common point of view though 
that the mode of energy injection in driven turbulence simula- 
tions is not important for statistical properties of the resulting 
flow dBoldvrev et al.ll2002h . The statistical analysis of our simu- 
lation data, however, cast doubt on his claim. Our results indicate 
that the dynamics of compressible turbulence at length scales 
significantly smaller than the energy-containing scale is aff'ected 
by the forcing. For this reason, the notion of fully developed 
turbulence-i. e., turbulence that is not affected by imposed con- 
straints such as boundary conditions, forcing or viscosity-is pos- 
sibly meaningless in the highly compressible regime. Of course, 
we cannot rule out that universal statistics might be asymptot- 
ically approached toward length scales smaller than the resolu- 
tion of our simulation. We note, however, that universality of the 
small-scale dynamics has been questioned even in the case of 
incom pressible turbulence dAlexakis et al.ll2005t ICasciola et al.l 
I2007h . 

The isothermal approximation is expected to apply to 
the interiors of molecular clouds, where the temperature of 
the cold ph ase of the interstellar ga s is close to its min- 
imum (Balles teros-Paredes et al.l [20071) . In molecular clouds, 
there are prominent magnetic fields. Although they have 
indisputably a large impac t on fragmentation properties 
jHennebelle & Tevssierl |2008|) . it appears that scaling expo- 
nents are not greatly altered by magnetic fields, at least in the 
super -Alfvenic case (Boldyrev et al. 2002t iPadoan et al.l 12001 
12007 ). Consequently, our analysis of the velocity statistics 
of compressively driven turbulence has potential relevance to 
molecular cloud turbulence. At larger scales in the ISM, var- 
ious heating and cooling processes bring about thermal insta- 
bilities. These are included in the colliding flow simula tions 
by [H eitsch et al. (2006), V azquez-Semadeni et al] (|2007|) and 
iHennebelle & Audit (.2007,) . We will consider driven turbulence 
in thermally bistable gas in a follow-up paper In the future, we 



also plan to include self-gravity and magnetic fields. By succes- 
sively including more complex physics, we expect to achieve 
a thorough understanding of the effects stemming from certain 
physical processes. The numerical treatment of processes such 
as cooling or self-gravity poses the problem of covering a wide 
range of length scales, at least in some regions of the flow. Using 
grid-based codes, a possible solution is adaptive mesh refine- 
ment (A MR). The simulation o f turbulence with AMR was pio- 
neered bv iKritsuk et al.l (l2006l) . Recently, an AMR simulations 
of self-gravitating turbulence have been presented (Offn eret al.l 
2008). As a preliminary study, we also performed an AMR sim- 
ulation, with the rate of compression (i. e., the rate of change 
of the negative divergence) and the vorticity modulus as control 
variables for refinement. Velocity and density statistics are well 
reproduced in comparison to the simulation with a static grid, 
but, using only one level of refinement, we have not achieved a 
significant gain in terms of performance in the case of isother- 
mal sugersonicturbu^^ The same difficulty was encountered 
by (I Kritsuk et al.ll2007ah for more than twice the resolution we 
used. 

The article is organised as follows. Section |2] explains nu- 
merical techniques, in particular, the stochastic forcing and 
the method of refinement employed in the AMR simulation. 
In Section [3] we begin our presentation of simulation results 
with global statistics and interpret some properties via three- 
dimensional renderings of the mass density and the vorticity in 
the stage of turbulence production. In the following, we anal- 
yse the probability density functions of the mass density and the 
vorticity. Next we turn our attention to the spectra of the tur- 
bulence energy and the density fluctuations. Finally, we derive 
scaling exponents from velocity structure functions. The results 
are related to each other and compared to observational results, 
theoretical predictions as well as results from other numerical 
simulations in Section H] In the conclusion, we discuss possible 
implications on molecular cloud turbulence and star formation 
and give an outlook to future work. 



2. Numerical techniques 

For the simulations pres ented in this article, we employed the 
open source code Enzo dO'Shea et al. I I2OO5I) in order to solve 
the compressible Euler equ ations by means of th e piec ewise- 
paraboHc method (PPM) of IColeUa & Woodward! (11984!) . The 
Euler equations including an external force density pf can be 
written as follows: 



D 
D 
D 

n — e + V ■ Pu - of ■ u. 



(1) 
(2) 
(3) 



The primitive variables are the mass density p, the velocity u 
and the specific total energy e of the fluid. The Lagrangian time 
derivative is defined by 



0-9 

= — -H M ■ V. 

Df dt 

The total energy per unit mass is given by 
1 



—u + 



2 {J -Dp 



(4) 



(5) 
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where y is the adiabatic exponent and the pressure P is related 
to the mass density and the temperature T via the ideal gas law: 



P = 



(6) 



The constants ^b, and mH denote the Boltzmann constant, the 
mean molecular weight and the mass of the hydrogen atom, re- 
spectively . 

The specific force / in equations Q and (O smoothly ac- 
celerates the fluid at large scales. Using a periodic domain, we 
adopt a generalisation of the Ornstein-Uhlenbeck process for the 
Fourier modes f{k, f) of th e acceleration field (lEswaran &Popel 
fT988t ISchmidt et al.ll2006h : 



1/2 



■ (7) 



In this stochastic differential equation, the Wiener process 'Wi 
generates Gaussian random vector deviates. We define the char- 
acteristic forcing wavenumber ko :- Ina/X for a cubic do- 
main of size X. Ideally, one would prefer to have a <K 1 a 
in order to minim ize the effect of periodic boundary conditions 
(lDavidsonll2004l) . However, this would constrain the dynamical 
range by far too much for a grid size that is computationally 
feasible. For this reason, we set ff = 2 so that the wavelength 
of the force field is about half of the domain size and there are 
~ 10 energy-containing structures ("large eddies"). To allow for 
a small spread of forcing wave numbers, we use the parabolic 
weighing function cr{k) oc k^(2ko - k)^H(k - Iko), where k - \k\ 
and H is the Heaviside step function, selects non-zero modes for 
1^1 €]0, 2A:o[. The autocorrelation time scale T of the forcing is 
identified with the large-eddy turn-over time, i. e., T - L/V, 
where V is the characteristic velocity of the flow in the station- 
ary regime and L - X/a - In/ko is the integral length. The 
stochastic diffusion term in equation O ensures that the result- 
ing force field becomes statistically isotropic in physical space, 
while any anisotropic initial condition decays exponentially by 
virtue of the first term on the right-hand side. This is not guaran- 
teed with the method of steady random forcing, where isotropy 
must be prepared to high precision in the initial velocity field. 

For purely solenoidal forcing, the random deviates generated 
by the Wiener process are projected perpendicular to the associ- 
ated wave vectors k by means of the projection operator 

(PuW) = iPfiik) + (1 - OPl(k) = ^6ij + (1 - 20^. (8) 

with ^ - 1. In the case ^ = 0, on the other hand, projection 
parallel to the wave vectors results in a dilatational forcing field. 
Due to the normalisation factor 



(9) 



in equation O, the root mean square acceleration f^ms = 3V^/L 
is independent of ^. We chose V - Scq/ V3, where cq is the ini- 
tial speed of sound. For the projection operator we set ( - 
0.1. Thus, the forcing in our simulations is constructed as super- 
position of a compressive component that produces supersonic 
flow with Mach number roughly given by Ma - V/cq ~ 2.9 and 
a subdominant solenoidal component. In consequence, there are 
two mechanisms of vorticity generation: Firstly, large-scale ed- 
dies are directly generated by the solenoidal forcing component. 



Secondly, the collision of shock fronts produced by the compres- 
sive component of the driving force resiilts in vorticity genera- 
tion due t o small-scale i nstabili ties (IWalder & Folinilll998l) . In 
contrast to lKritsuk et al.l (l2007ah . the second mechanism is more 
important in our simulations. 

In the limiting case that all heat produced by kinetic en- 
ergy dissipation is removed from the system immediately, the 
gas evolves isothermally. Nearly isothermal turbulence can be 
simulated by setting the adiabatic exponent y to a value slightly 
greater than unity. Then P ^ pc^, where Cj is approximately 
equal to the constant isothermal speed of sound. For j - 1.01 
and Ma^ < 10, as in our simulation, the specific kinetic energy 
Ckin ^ lOCj and the internal energy density ejnt - P/ij - ~ 
lOOCj . Thus, the dissipation of kinetic energy results in a secu- 
lar change of the internal energy. In this sense, the g as is nearly 
isothermal. As pointed out by Kritsuk et al. (2007a), the mean 
helicity of isothermal turbulence should be conserved. In gen- 
eral, the velocity increments produced by random forces are not 
free of helicity locally but average out globally. The helicity of a 
dilatational force, however, is identically zero. 

The main simulation run on a static grid of N = 768^ 
cells, while a root grid of A^o - 196^ and one level of refine- 
ment with an increase of resolution by a factor 4 was applied 
in an AMR simulation. The effective resolution, A^i,eft = 768^, 
matches the resolution of the static grid run. An essential re- 
quirement for computing turbulent flow with AMR is the sensi- 
tivity of the refinement criteria on small-scale properties of the 
flow. For this reason, gradients of the velocity, the density or 
other fields are commonly used as control variables for refine- 
ment (Kritsuk et aL 2007a). In this article, we focus on struc- 
tural invariants, i. e., scalars derived from the velocity gradi- 
ent. A choice that suggests itself for tracking down turbulent 
eddies is the enstrophy. Thus, we define the control variable 
hi - ju/, where the vorticity cj - V x m. In supersonic tur- 
bulence, the steepening of density gradients is associated with 
rising gas compression. This is indicated by a positive rate of 
convergence, i. e., the negative time derivative of the divergence 
d - V ■ u. Letting the operator V- act upon equation (|2]), it can 
be shown that the rate of convergence is given by 

-§-d^l{\Sf-^^).V.[lvP-f]. (10) 

In the isothermal limit, this equation becomes 

-^^d=^-{\S\^-aj^) + clV^6-V-f (11) 

where 6 = ln(p/po)- Dropping the large-scale forcing term on 
the right hand side, which is not supposed to trigger refinement, 
we define the control variable 



(12) 



Note that /!2 - in the incompressible limit. In more gen- 
eral applications such as turbulence in non-isothermal or self- 
gravitating gas, the corresponding effects can be singled out as 
additional terms in /j2. 

Regardless of the refinement variables in use, one has to 
find an appropriate normalisatio n in order to set thr esholds for 
triggering refinement. Whereas Kritsu k et al. I (l2007a>) normalise 
the control variables in their AMR simulation with the local 
flow velocity and gas density, we determine thresholds for re- 
finement from statistical properties of the associated control 
variables. This has the advantage that no a priori knowledge 
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Table 1. Global statistics at various instants of time. From left to right, the columns contain the normalised time, peak density, 
standard deviation of the density, cumulative probability of Mach numbers > 1, RMS Mach number, RMS vorticity, mass- weighted 
RMS vorticity, Taylor scale, mass-weighted, Taylor scale, small-scale compressive ratio and the corresponding mass-weighted 
compressive ratio 



about the control variables is required. For the calculation of 
the thresholds, both the spatial average (/?,(jc, t)) of the control 
variable h-, and its standard deviation std hi{t), where std^/i,(f) — 
{h^{x, 0) - {hiix, tyf-, are used. A cell will be flagged for refine- 
ment, if the fluctuation 

h'. = hi - {hiix, t)) > thresh /!,(f), (13) 

where the threshold is defined by 

thresh /z,(f) := A max [{hiix, f)>, std /i,(f)] ■ (14) 

This definition of the threshold ensures that refinement is only 
triggered by pronounced fluctuations. For a nearly uniform 
field with h'-ix, t) <K {hiix, t)), the above criterion amounts to 
h'-ix, t) > A{hiix, t))i, i. e., virtually no cells are refined provided 
that A > 1 . For the opposing limit of a field with large fluctu- 
ations relative to the average, the refinement criterion becomes 
h\ix, f) > A std hiit). In the general case of multiple levels of re- 
finement, the averaging has to be constrained to refined regions 
or grid patches at higher levels. The efficiency of refinement can 
be tuned by the coefficient A . If A » 1 , only peak fluctuation will 
trigger refinement. From test simulations with various choices of 
A, we concluded that statistical properties of turbulence are well 
reproduced for A = 1, as one would intuitively expect. This set- 
ting that was used for ffie production run. 

3. Simulation results 

Making use of the invariance properties of isothermal turbu- 
lence, we normalise all variables in the following such that the 
mean density po = 1, the initial pressure Pq = 1 and the size 
of the computational domain X - 2L - 1. This corresponds 
to the internal representation of the fluid-dynamical variables in 
Enzo. The speed of sound is ffien given by cq = yjy ~ 1.005, 
the characteristic velocity V ^ 2.901 and the integral time scale 
T X 0.172. The static-grid simulation run for < f < 1.6, span- 
ning about 9 integral times scales. Output was produced in in- 
tervals of 0.025 resulting in 7 data dumps per integral time. An 
overview of global mean quantities is given in Table [T] For an 
illustration of the evolution within the first integral time scale, 
ffiree-dimensional volume renderings of the mass density p and 
ffie vorticity modulus o) at time tjT - 0.44, 0.58, 0.73 and 1.02 
are shown in Figures [T] and |2l respectively. These volume ren- 
derings were produced from AMR data. The AMR simulation 



was stopped when the first integral time scale was completed be- 
cause, at this point, it progressed slower than the static-grid run. 
In part, this was caused by a severe performance drop of Enzo as 
the total number of grid patches approached lOOOOQ Apart from 
that, more than 60 % of the computational domain is covered by 
refined grid patches at time tjT ^ 1 . With a single level of refine- 
ment, it is difficult to compete with a static-grid simulation for 
such a volume filling factor The development of the mesh struc- 
ture in two-dimensional slices is illustrated in Figures [3] and |4] 
Later in this Section, we will show that statistical properties in- 
ferred from both runs agree very well. In the conclusion of this 
article, we will comment on perspectives of applying AMR to 
the numerical simulation of turbulence. As regards the evolution 
over several integral time scales, we will refer to the static-grid 
data in the following. 

Initially, as one can see from the maximal values and the 
standard deviations in Table[Tl the mass density rises rapidly and 
reaches a peak roughly at half an integral time scale. The RMS 
Mach number, Ali-ms = ((v/cs)^)'^^, is then close to 3. The reason 
for the fast production of density enhancements of order 10^ is 
illustrated by panels (a) and (b) in Figure[T] The outflows of gas 
from proximate regions, in which the divergence of the stochas- 
tic force is positive, collide and become supersonic on a time 
scale ~ r/Ma ^ Q.2T. Interpolation of the cumulative probabili- 
ties for Mach numbers greater than unity implies that about 60 % 
of the domain is filled by supersonic flows at time t - 0.2T. The 
collision of these supersonic flows results in the sheetlike struc- 
tures of shock-compressed gas that can be seen in the volume 
renderings. The density structure is compared to the Mach num- 
ber Ai = v/cs of the gas in the two-dimensional slices plotted 
in Figures |3] and |4] Since the forcing evolves on a slower time 
scale than the colliding flows, gas accumulates in the over-dense 
regions faster than the drift of the positive-divergence regions 
could act against further compression. On the other hand, com- 
pressed gas re-expands in between reflected shock fronts travel- 
ling in opposite directions. Presently, the gas becomes turbulent 
due to instabilities at the collision interfaces and the solenoidal 
forcing component. As one can see in Figures [3] and [H these in- 
terfaces are well tracked by AMR, while smooth regions of over- 
dense gas are not generally refined. At the end of the first inte- 
gral time, cloud-like structures have emanated from the vertices 
of the high-density sheets (see panels (c) and (d) in Figures [TJ. 

' Meanwhile, this problem has been resolved by the Enzo developer. 
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(a) t = 0.44r, m = 429&,xt = 0.299 



(b) t = 0.5ST, m = 5208, ;ri = 0.321 




(c) t = 0.73r, «i = 6331,^1 = 0.423 



(d) t = 1.02T, n, = mi5,xi = 0.653 



Fig. 1. Volume renderings of the mass density p for different stages of developing turbulence in the AMR simulation. For each 
instant of time, the total number of grid patches at the first level of refinement, A^i, and the corresponding volume filling factor 
are specified. 



The maximal densities are somewhat smaller at this stage and 
A^i-ins ~ 2.5. Subsequently, A^i-ms slowly decreases because of 
the gradual increase of the internal energy caused by energy dis- 
sipation. 

A crucial question for the following analysis is whether tur- 
bulence settles into a steady state, in which small-scale scale ve- 
locity fluctuations are fully developed and kinetic energy dis- 
sipation is balanced by the integral-scale supply of energy from 
the forcing. To this end, we consider several statistical measures. 



The Taylor scale A follows from the two-point autocorrelation of 
the turbulent velocity field. For incompressible isotropic turbu- 
lence, A is given by the the time scale associated with th e RMS 
vortic ity, l/oj^ms, times the RMS velocity fluctuation dFrischl 
[1991 : 

^=V5 — . (15) 

^nns 

Using the above expression as definition of A also in the case 
of compressible turbulence, we calculated the values listed in 
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(a) t = 0.44r, m = 429&,xt = 0.299 




(c) t = 0.73r, «i = 6331,^1 = 0.423 




(b) t = 0.5ST, m = 5208,^-1 = 0.321 




(d) r = 1.02r, n, = mi5,xt = 0.653 



Fig. 2. Volume renderings of the vorticity modulus oj corresponding to Figure[T] 



Table [T] Initially, A ~ L, indicating that the flow is smooth and 
only varies over the integral length L. In the course of the first 
integral time, the Taylor scale drops substantially and then grad- 
ually approaches the asymptotic value A x 20. The early evo- 
lution of A reflects the growth of small-scale structure that can 
be seen in Figure |2] As velocity fluctuations develop at smaller 
scales, the vorticity rises rapidly and, hence, A decreases. 



Based on the equilibrium between energy injection and dis- 
sipation in the statistically stationary state, the Taylor scale can 



(16) 



be related to the Reynolds number (|Popell2000h : 

A _ fm 

Z ~ V Re' 

This relation ofifers the possibility of estimating the Reynolds 
number from statistical properties of the velocity field without 
explicit knowledge of the numerical viscosity. Substituting L = 
384A and A/ A 20, we obtain Re a; 3.7 ■ 10^ Since the criterion 
Re = {LjKfl^ yields Re » 2.8 ■ 10^, the Reynolds number in our 
simulation is comparable to the Reynolds number that could be 
achieved in a simulation of incompressible turbulence with the 
same dynamical range. 
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(c) p(t = 0.58r) (d) M(t = 0.5ST) 

Fig. 3. Contour plots of the mass density p and the Mach number Ai in the plane z = corresponding to the top panels in Figures [T] 
and|2] The rectangles show the boundaries of refined grid patches. 



Since the Taylor scale only probes the vorticity of the veloc- 
ity field, we also consider the so-called small-scale compressive 
ratio (.Kida&Orszag..l990.) . 



(17) 



At early time, this ratio is close to unity, i. e., the rotation-free 
component of the velocity dominates. This reflects the large- 
scale stochastic force. Subsequently, however, vorticity develops 
in the vicinity of shock-fronts and grows due to the production of 



turbulence. The compressive ratio in the steady state is slightly 
less than 0.5. Thus, the solenoidal and dilatational modes of the 
velo city fluctuations are co mparable in magnitude. For compar- 
ison, |&itiuketal] (|20073) found ~ 0.28 in their simulation. 



The definitions of the Taylor scale ( fTSl l and the small-scale 
compressive ratio ( fT9] l stem from the theory of incompressible 
turbulence. We also calculated mass-weighted variants based on 
the kinetic energy density, jpu^, the density of the enstrophy. 
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(c) pit = 1.02T) (d) M(t = 1.027) 

Fig. 4. Contour plots of the mass density p and the Mach number Ai in the plane z = corresponding to the bottom panels in 
Figures [1] and |2] The rectangles show the boundaries of refined grid patches 



and jpd . The values of the mass-weighted Taylor scale. 



(18) 



listed in Table [T] show a faster decline in the course of the first 
integral time in comparison to A. This suggests that a significant 
fraction of the vorticity originates from shock-compressed gas. 
Indeed, panels (a) to (c) in Figure |2] show sheet-like structures 
which can be attributed to shock fronts at the interface between 
compressed and rarefied gas. Moreover, Figure |4] suggests that 



the production of turbulence ensues in the vicinity of shocks. For 
tjT > 2, we have A^^jA x 0.6. Note that the intermittency of 
the mass density causes stronger fluctuations of Amw- The mass- 
weighted small-scale compressive ratio. 



(19) 



asymptotically approaches a value slightly less than 0.3 « 0.6rcs, 
indicating that solenoidal flow is somewhat more pronounced 
relative to dilatational flow in over-dense gas. Altogether, the 
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Ttj 



Fig. 5. Temporal evolution of the probability density functions 
of the vorticity modulus oj. 




InO/po) 

Fig. 6. Instantaneous and temporally averaged (thick solid line) 
logarithmic mass density fluctuations, ln(p/po). 



evolution of the Taylor scale and the small-scale compressive 
ratio shows that it takes about two time scales for the turbulent 
flow to become statistically stationary, possibly, with small sys- 
tematic changes up to t/T ^ 4. 

In the remainder of this Section, the statistics of the mass 
density and the turbulent velocity field will be analysed further 
by means of probability density functions, spectrum functions 
and structure functions. 

3. 1 . Probability density functions 

For statistical analysis, turbulent flow quantities such as the gas 
density can be considered as spatiotemporal random variables. In 
numerical simulations based on finite-volume schemes, Eulerian 
statistics is most readily obtained by calculating probability dis- 
tributions for the whole computational domain at any given in- 



stant of time. The probability distribution function (PDF) quanti- 
fies the cumulative probability that a random variable is found to 
be less than a certain value. The derivative of the probability dis- 
tribution function with respect to value of the random variable 
is called the probability density function (pdf)Q In the follow- 
ing, velocity and mass density statistics is specified by means of 
probability density functions. 

The pdfs of the vorticity modulus plotted in panel (b) of 
Figure|5]highlight the rapid evolution of the flow within the first 
integral time. At f » 0.5T, a power-law decrease of the w-pdf 
in the subrange 10 < w < 100 becomes apparent and flattens 
subsequently. The power-law behaviour points at a self-similar 
ensemble of shocks. The flattening arises form the formation of 
small-scale structure, as illustrated in Figure|2] because velocity 
fluctuations at smaller scales are associated with higher vortic- 
ity. At later time, the power-law subrange diminishes and the 
exponential tail towards high vorticity, which is characteristic 
for incompressible turbulence, becomes more prominent. This 
indicates a transition from an early shock-dominated phase into 
vortex-dominated turbulence in the course of the first few inte- 
gral times. From t ^ 3T onwards, the w-pdfs show only little 
variation in time. 

The logarithmic density fluctuations, 6 - ln(p/po), on the 
other hand, show no systematic evolution in the aftermath of 
the shock-dominated production phase. This corresponds to the 
nearly constant RMS Mach number for t/T > 1 (see Table [TJ. 
In Figure |6] (5-pdfs are plotted for several instants of time to- 
gether with a time-averaged pdf. For the averaging 57 data sets 
produced in the time interval 1 < t/T < 9 were used. One can 
see random variations of the instantaneous pdfs, which reflect 
intermittent deviations from the ensemble average. The average 
of the pdfs over several integral times, on the other hand, yields 
a sound approximation to the statistics of the mass density in 
the ensemble average. Most noticeably, the averaged (5-pdf is 
not symmetric but exhibits negative skewness, i. e., lower densi- 
ties are more probable in relation to higher densities. This corre- 
sponds to the appearance of compact regions of over-dense gas 
surrounded by extended voids in Figure[T] 

Probability density functions are indicators for the reliabil- 
ity of AMR simulations, because pdfs would be affected, if the 
fluctuations of the corresponding quantities were not properly 
tracked by refinement. Here, we consider pdfs of the mass den- 
sity, the vorticity modulus and divergence. For the termination 
time of the AMR simulation (f - 1.027), Figures |8] [T] and 
|9] show that the pdfs calculated from the AMR data are well 
matched by the corresponding static-grid pdfsQin particular, the 
high-density and high-vorticity tails arising, respectively, from 
strongly compressed gas and from thin vortex filaments are very 
well reproduced with AMR. There are deviations in the left tails 
of the pdfs though. This is to be expected because the refinement 
criteria are optimised to capture peak values that contribute to 
the tails on the right. The important physics is associated with 
these tails. Since turbulence is close to statistical equilibrium af- 
ter the first integral time scale, the comparison demonstrates that, 
in principle, the AMR techniques used in our simulation are ap- 
plicable to turbulence simulations. For conclusive testing, how- 
ever, an improved code version is necessary to carry out AMR 
simulations with multiple levels of refinement over several inte- 
gral time scales. 



- Here, we adhere to the convention that the symbol of the derivative 
is written in lower case letters, while the primitive function is written in 
upper case letters. 

^ Note that the pdf of p is related to the pdf of S via p pdf (p) = pdf ((5). 
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Fig. 7. Comparison of the probability density functions of the 
vorticity oj for the static grid and the AMR run at time t - I .Q2T. 




Fig. 9. Comparison of the probability density functions of the 
divergence d for the static grid and the AMR run at time t = 
I.Q2T. 




0.001 0.010 0.100 1.000 10.000 100.000 
p/Po 

Fig. 8. Comparison of the probability density functions of the 
mass density p for the static grid and the AMR run at time t - 
1.02T. 



3.2. Turbulence energy spectra 

The turbulence energy spectrum function (lFrischlll995h . 

E{k) = ^ ^dQi k^u(k) ■ u{k) (20) 

where u{k) is the Fourier transform of the velocity field and u*{k) 
its complex conjugate, is the surface integral of the spectral en- 
ergy density per unit mass over a sphere of radius k in Fourier 
space. In the ensemble average, the energy spectrum function 
of statistically stationary turbulence is independent of time and 
obeys a power law, E{k) oc k^P, for fco - 2n/L <K A: <sc A:c, where 
kc = 7t/A is the numerical cutoff wave number We computed 
discrete energy spectrum f unctions form several data sets fol- 
lowing the prescription of iSchmidt et al.l (l2006i) . The resulting 
spectrum functions were compensated with the factor k^ so that 



the compensated spectrum of Burgers turbulence appears con- 
stant in the inertial subrange. Moreover, the normalization k - 1 
for the smallest wave number (corresponding to the box size) is 
applied. 

The sequence of compensated energy spectrum functions for 
consecutive instants of time shown in Figures [TO] and (TT] illus- 
trates the gradual build-up of the turbulent energy cascade in 
the course of the first integral time scale followed by the re- 
laxation of the flow into a steady state. In the following, it is 
understood that the spectral indices /3 are obtained by fitting 
power laws k/^ to the compensated spectrum functions, where 
/?' - p - 2. The spectra for t ^ T are substantially steeper than 
the Kolmogorov spectrum E{k) oc k^^^^. In particular, we find 
/? ^ -2.02 + 0.06 at time f = 1 .02 T, which is very close to /? = 2 
for Burgers turbulence and verifies an early shock-dominated 
phase. Subsequently, the spectrum functions become flatter and 
a bump roughly centred at a; 70 ~ fcc /5 appears . This bump, 
which stems from th e so-called bottleneck effect dPobler et al.l 
!20Q3t ISchmidt et al.l 1200 6), severely constraints the power-law 
range of the spectra. Applying, least square fits in the wave num- 
ber range 5 < A: < 15 yields the time-averaged estimate /? a; 1 .87 
for f/r > 2 (see Table|2]i. However, due to the gradual flattening 
toward higher wave numbers caused by the bottleneck effect and, 
possibly, other factors related to PPM (iKritsuk et al.ll2007al) . the 
systematic uncertainty of this estimate is much higher than the 
fitting error margins. 

The integral of E{k) over the range of numerically resolved 
wavenumbers is the mean kinetic energy per unit mass; 

P E{k)dk = (21) 

where a - X/L and ko - 2nlL (see Section |2]i. For supersonic 
isothermal turbulence, Wnns > cq. Hence, there exists a wave 
number k^ > ko/ a, for which the cumulative turbulence energy 
in the range k > ks equals j{c^) - Cq, i. e., 

J ' Eik)dk - i<c,)2. (22) 
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Fig. 10. Early evolution of the energy spectra for t < T. The spectrum functions are compensated with k^. The sonic wave number 
is indicated by vertical dotted lines. Additionally, the Kolmogorov slope -5/3 is depicted as dashed line. 



The sonic wave number A;, corresponds to the length scale (g — 
In/ks, for which the turbulent velocity fluctuations m'(4) ~ Cs- 
The values of obtained by the inversion of equation ( l22b are 
marked by vertical lines in the plots of the energy spectrum func- 
tions (Figures [TOlandfTTI). In the course of the first integral time, 
kg reaches a maximum of about 3kQ. Together with the maximum 
of kg comes the steepest slope with a spectral index of /3 ~ 2.0. 
Afterwards, the sonic wave number settles at 2ko. One should 
note that no implications can be made with regard to the scal- 
ing for k > kg. Although it would seem that this subrange of 
wave numbers is dominated by nearly incompressible modes, 
the scaling is probably not consistent with the Kolmogorov the- 



ory. Since the uncertainty in the determination of /3 is quite high 
due to the bottleneck effect, a stronger case for this conclusion 
will be made on grounds of the exponent of the second-order 
structure functions. 

3.3. Mass density spectra 

Analogous to O(^) = u{k) ■ u{k)*, which specifies the density of 
specific kinetic energy, ^u^, in Fourier space, the spectral density 
of the squared logarithmic density fluctuation 6^ is given by 

<:>e(k) = m ■ 6*(k), (23) 
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Fig. 11. Temporal evolution of the compensated energy spectra continuing FigurefTOlfor t>T. 



and, by virtue of Parseval's theorem. 



(24) 



For statistically isotropic turbulence, the spectral density Q> f(k) 
of some field variable / is a function of the wave number k - \k\ 
only. The spectral density O^(fc) is the counterpart to <i>{k) in as 
much as both are related to squared fluctuations of the respective 
field variable. The numerically computed spectral density <i>s{k) 
compensated with k* is plotted in Figure[T2] Similar to the turbu- 
lence energy spectra, there is only a marginal inertial subrange. 
Form least-square fits of the compensated spectral density ^s{k) 
in the range 10 < A; < 30, we obtained the spectral indices /3g 



of the quadratic logarithmic density fluctuations. Averaging the 
indices listed in Table|2]for t/T x 2 to 6 yields j3s ~ 1.80. 



3.4. Velocity structure functions 

The turbulence energy spectrum function E{k) specifies the scal- 
ing behaviour of the squared velocity fluctuation in Fourier 
space. The corresponding description in physical space is given 
by the second-order structure function, i. e., the spatial two-point 
correlation function of the velocity field. Structure functions of 
order p are defined by 



Spi{) = {\u{x + {)-u{x)n ~{^", 



(25) 
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t/T 


p Pli 


0.73 


1.99 ± 0.10 


0.87 


2.04 ± 0.05 


1.02 


O AO 1 A A/; 

/.Uz ± O.Oo 


Z.Ui 


1 OA 1 A AC 1 OC 1 A A/; 

1.89 ±0.05 1.83 + U.Uo 


3.05 


1.85 ±0.05 1.76 ±0.04 


4.06 


1.92 ±0.04 1.75 ±0.05 


4.93 


1.87 ±0.08 1.82 ±0.05 


5.95 


1.90 ±0.04 1.80 ±0.04 



Table 2. Spectral indices of turbulence energy and mass density 
fluctuation spectra with fit error margins. 
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Fig. 12. Compensated spectral density of the quadratic logarith- 
mic density fluctuation at time t/T - 5.95. 



where i - In/k specifies the length scale. Of particular im- 
portance are the exponents ^2 and ^3 which are related to the 
scaling of turbulence energy and the turbulence energy flux, re- 
spectively. |KritsuketalJ(l2007a|) proposed to calculate structure 
functions for the mass-weighted velocity p'^^H. This is moti- 
vated by the dependence of the energy flux on pu^ in compress- 
ible turbulence. However, observational results for the scaling 
of molecular cloud turbulence are only available from measure- 
ments of the velocity dispersion at different scales. For this rea- 
son, we restrict our discussion to the scaling properties of the 
structure functions without mass-weighing dZSl l. 

The computation of statistically converged structure func- 
tions of higher order is a computationally demanding task. We 
applied a Monte-Carlo algorithm to sum velocity differences 
from randomly sampled points over bins of spatial separation. 
All results presented hereafter have been computed based on ap- 
proximately 2.4 X 10'*' data pairs. To ensure that sufficient num- 
bers of data pairs were sampled, we increased the amount of 
sampled data by a factor of 5 in one particular case. No sig- 
nificant differences were found for the structure functions up to 
fifth order. Since the measurement of line broadening at posi- 
tions separated by a certain angular distance within molecular 
clouds corresponds to transversal velocity fluctuations, we com- 
puted transversal structure functions, 5y,,tians(^), for which the 
velocity increments u(x + f)-u are projected perpendicular to {. 
For t/T a: 1, 2, 3 and 4, the structure functions of second order 
and third order are plotted in Figure [T3] We fitted power laws 
in the subrange 8 < ^/A < 50 for t/T 1 and 8 < (/A < 70 
for t/T > 2, respectively, which are also shown in Figure [T3] 



Table [3] lists the resulting scaling exponents for p = 1, . . . , 5 
in intervals of about one integral time scale. As one can see, 
52,trans IS Steeper at time t/T a: 1 in comparison to later instants 
of time. The scaling exponent (2 ~ 0.95 is close to the exponent 
of Burgers turbulence. For t/T > 2, there appears to be little 
variation of 52,trans, although there might be a slight systematic 
decrease of the scaling exponents till t/T a 4. 

Plotting 52,trans against 53_trans, as shown in panel (a) of 
Figure [141 reveals power-law behaviour that partially extends 
into the dissipation range and also applies to the largest scales 
i ~ L. This remarkable property is known as extended self- 
similarity (ESS) (.Benzi et al 1993.). For t/T a 2, the struc- 
ture functions 5p,trans are shown as functions of ^a.trans for 
p - 1, . . . , 5 in panel (b) of Figure [14] It appears that ESS is 
very well satisfied for the transversal structure functions with 
p < 5, albeit with a somewhat smaller range in the case p = 5. 
The slopes of the corresponding power-law fits functions for 
0.04 < S'B.trans < 1 -0 yield relative scaling exponents Zp = (p/(-i. 
Averaging the values of Zp for t/T > 2, we obtained the values 
that are summarised in Table |4] The averaged scaling exponents 
are plotted with the corresponding scattering range for each p in 
Figure [Is] 



4. Discussion 

4. 1 . Velocity statistics 

The second order exponent ^2 determines the inertial-range scal- 
ing of the RMS velocity fluctuations, u'{{) oc (^^ where y - (2/2. 
From the power-law fits to transversal structure functions of 
second order for t/T > 2, we obtain the average ^2 ~ 0.868 
corresponding to y a 0.43, while y a 0.48 is implied for 
t/T X 1. The spectral index - 2j + I (see Section [J!2] i im- 
plies y ss 0.5 in the transient phase prior to equilibrium and 
y Si 0.45 for statistically stationary turbulence in our simula- 
tion. Given the large uncertainty of fitting the turbulence en- 
ergy spectra, the spectral indices are consistent with the scal- 
ing exponents obtained from the second-order structure func- 
tions. Observations of molecular clouds indicate power law re- 
lations between the velocity dispersion and the size, Av oc 
The index a is interpreted to correspon d more or less to y, al - 
though there is not a s trict equivalence (iBrunt & Heveiil2002l) . 
While Falgar one et al.l (.1992) found a ^ 0.4 and a ^ 0.43 
was obtained by Miesch & Ballvl (1 19941) in very good agree- 
ment with the steady-state v alue of y inferred from our simu- 
lation, iBrunt & Heverl (l2002h deduced from their measurements 
a much higher mean y a 0.57. However, their results for indi- 
vidual molecular clou ds vary in the range y a 0.33 . . .0.81. For 
Perseus. [Padoan et al.l ([2006) obtained an index /3 ~ 1.81 corre- 
sponding to y a; 0.4, again very close to the index we estimated 
for the turbulence energy spectra in statistical equilibrium. 

Other than the numerical determination of spectral indices 
and absolute scaling exponents, relative scaling properties can 
be inferred quite reliably (see Figure [l4] l. In Table [4] various 
theoretical predictions and numerical results for the relative scal- 
ing exponents Zp = (p/^i of order p - 1, . . . , 5 are listed. The 
classical Kolmogorov theory of statistically stationa ry isotropic 
incompressible turb u lence yields Z^ - ^"p = p/3 (.Kolmogorov! 



mcompressible turp u lenc 

119411; lFrischlll995h . IShe & Levequg (1 19941) included intermit- 
tency corrections for incompressible turbulence based on the 
idea of a hierarchy of dissipative structures. The scaling expo- 
nents predicted by t heir model were experimentally confirmed 
with high precision. lDubrullel ( ll994 showed that this model is 
a member of the general class of hierarchical structure models 
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tIT 


P = 


1 


P 


= 2 


P = 


3 


P = 


4 


P 


= 5 


1.02 


0.654 ± 


0.002 


0.950 


± 0.003 


1.107 ± 


0.006 


1.202 ± 


0.009 


1.256 


±0.012 


2.03 


0.562 ± 


0.004 


0.888 


± 0.003 


1.066 ± 


0.008 


1.167 ± 


0.013 


1.228 


±0.016 


3.05 


0.538 ± 


0.005 


0.884 


± 0.003 


1.097 ± 


0.006 


1.239 ± 


0.013 


1.344 


±0.019 


4.06 


0.531 ± 


0.003 


0.858 


± 0.004 


1.035 ± 


0.013 


1.123 ± 


0.023 


1.163 


±0.031 


4.93 


0.530 + 


0.002 


0.848 


± 0.006 


1.018 ± 


0.016 


1.100 ± 


0.026 


1.135 


± 0.034 


5.95 


0.523 + 


0.003 


0.852 


± 0.004 


1.044 ± 


0.012 


1.157 ± 


0.021 


1.224 


± 0.027 


6.96 


0.546 + 


0.004 


0.875 


± 0.004 


1.053 ± 


0.010 


1.150± 


0.017 


1.208 


± 0.021 



Table 3. Scaling exponents obtained from power-law fits to the structure functions 5p,trans(0 for several stages of the turbulent 
flow evolution. 



based on log-Poisson statistics. Consider the local rate of dissi- 
pation £( at some length scale {. Going to smaller length scales, 
the rate of dissipation tends to grow as structure at smaller scales 
emerges. This implies the scaling law €( oc where A is the 
scaling exponent of the most singular dissipative structures as- 
sociated with the divergent rate of dissipation as f — > oo. On the 
other hand, intermittency tends to decrease the rate of dissipa- 
tion. Let y6 be a parameter (not to be confused with the spectral 
index of turbulence energy) that speci fies the degre e of inter- 
mittency in the fashion of the /3 model ('F risch|[T995l) and let us 
assume that the dissipation rates at different length scales are 
connected by n events out of a Poisson distribution, which break 
up the dissipative structures from larger toward smaller length 
scales. Then, for any two length scales £2 < (i. 



(26) 



Because of the law of finite energy dissipation dFrischl [19951) . 
we have the normalisation constraint (e^-,) = (e^,) = e, where 
e is the mean rate of energy dissipation. Invoking the Poisson 
distribution, it follows that (e^^_oc -^j^j^ ^j^g 

refined similarity hypothesis (lFrischlll995h . the relative scalings 



Z, = (l-A)f.^(l-/.W3) 



(27) 



for th e velocity structu re functions ar e obta ined (iDubruUd 
[1991)0 To determine A, IShe & Levequd (11994 postulated a a 



uniform time scale oc for the dissipation of various turbulence 
intensities. The scaling of the most singular dissipative structures 
is given by the largest available kinetic energy, which is of the or- 
der y^, divided by l'^. In the case of incompressible turbulence, 
A = 2/3. Setting /? = 2/3, they interpreted C = A/(l -/?) = 2 
as codimension of one-dimensional vortex filaments, which are 
considered to be the most singular dissipative structures of in- 
compressible turbulence. 

For supersonic turbulen ce, in which dis sipation is attributed 
to two-dimensional shocks, BoldvrevI (l2002 l) proposed to set C - 
1 while leaving the scaling index A — 2/3 unchanged. Then 
P - 111) and, hence. 



Z„ = ^ + 1 - 



(28) 



The predictions of this so-called Kolmogorov-Burgers model for 
Z(l) and Z2 agree very well with numerical results obtained 
from magnetohydrodynamic turbulence simulations with purely 
solenoidal forcing ( Boldvrev et al. 2002) an d the supersonic tur- 
bulence simulation bv lKritsuk et alj (l2007al) . However, the time- 



Our simplified motivation follows [Panet '^ (I2008h . 



averaged values of Zp obtained from our simulation differ sig- 
nificantly the Kolmogorov-Burgers model for all p (see Table |4| 
and Figure [TSl ). IKritsuk et al.l (j2007b) found deviations of the 
scaling exponents of mass-weighted velocity structure functions 
for p > 3. 

There are several factors, from which the discrepant scaling 
exponents might stem. Firstly, it is possible that the computation 
of the structure functions is not well converged due to insuffi- 
cient sampling. Especially, this might affect the scaling expo- 
nents of higher order. However, varying the sample size did not 



W. Schmidt et al.: Compressively driven turbulence in isothermal gas 



15 




Fig. 14. Transversal structure functions plotted against the corresponding third-order structure functions. Panel (a) shows 52,trans as 
function of ^s uans for the same instants of time as in Figure [T3] In panel (b), structure functions (thick lines) and the corresponding 
power-law fits (thin lines) up to fifth order at time t » 2.Q3T are shown. 



p 


K41 


SL94 


B02 


BNP02 


KNPW07 


this work 


D94 


D94 






(C = 2) 


(C= 1) 


(MHD 500^) 


(HD 1024^) 


(HD 768') 


(C=0.71) 


(C=1.23) 


1 


0.33 


0.36 


0.42 


0.42 


0.43 


0.52 


0.54 


0.53 


2 


0.67 


0.70 


0.74 


0.74 


0.76 


0.83 


0.82 


0.83 


3 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


4 


1.33 


1.28 


1.21 


1.20 




1.09 


1.14 


1.10 


5 


1.67 


1.54 


1.40 


1.38 




1.14 


1.26 


1.16 



Table 4. Comparison of the normalised scaling exponents Z p = for structure functions of order p — 1 , . . . , 5 predict ed by 

iKolmogorovl (|1941|) . IShe & Levequd (11994 . iBoldvrevI (|2002|) and the general hie r archical structure mod el of lDubrullj (11994 ') with 
numerical results for transversal structure functions from Boldvrev et al.l (|2002|) . iKritsuk et al.l (l2007al) and this work. Note that 
Mms ~ 10 for BNP02 and Mr^s ~ 6 for KNPW07. 



affect the obtained scaling exponents noticeably, so we consider 
the calculated structure functions to be converged. Secondly, the 
resolution of our simulation is possibly not high enough to al- 
low for reliable estimates of the scaling exponents. Other than 
the turbulence energy spectra, however, structure functions show 
a relatively broad power-law range since there is no bottleneck 
effect (see Table[3]and Figure [T3li and ESS can be utilised to cal- 
culate the relative scaling exponents even more reliably. Already 
from 256-' test runs, we were able to obtain estimates of the rela- 
tive scalings that showed the same trend as in the high-resolution 
simulation. Thirdly, we consider scaling exponents of transver- 
sal structure functions, while theoretical model usually refer 
to longitudinal structure functions. However, we found similar 
values of Zi and Z2 from the longitudinal structure functions, 
while the disagreement of the higher-order scaling exponents 
with the Kolmogorov-Burgers model appears to be even more 
pronounced. Fourthly, the statistics might not be converged in 
time. This is excluded on grounds of the global statistics anal- 
ysed at the beginning of Section[3]and the small variations of the 
instantaneous scaling exponents listed in Table [3] There is pos- 
sibly a small systematic drift of the early scaling exponents, but 
this has no significant influence on the results discussed in the 
following. 

Given the relatively low RMS Mach number (A^rms - 
2.2... 2. 5) in our simulation, it is possible that vortices con- 
tribute more in comparison to shocks compared to the other 
simulations cited above, for which Alims > 5. On grounds of 
results from MHD simulations with different Mach numbers. 



i Padoan et al.l (|2004|) proposed a one-parameter model based on 
equation ( |27] i, where j3 = I - 2/(3C) and the the codimension 
C should monotonically decrease with Alims from C = 2 in the 
subsonic limit to C = 1 for sufficiently high RMS Mach num- 
bers: 



1 - 



2 

3C 



(29) 



The best fit of this equation to the time-averaged scalings in our 
simulation yields C ~ 0.71, which falls outside the range of the 
one-parameter model. However, one should keep in mind that 
the corresponding fractal dimension D » 2.3 might very well 
result from the stretching and folding of dissipative structures in 
turbulent flow. Letting p vary continuously, the corresponding 
curve is plotted in Figure [15] together with the data points and 
other model predictions. 

If the notion of universal scaling in the sense of the one- 
parameter model was correct, then, particularly, the second- 
order exponent Z2 would vary monotonically with the RMS 
Mach number form 0.7 in the subsonic limit to 1 .0 in the hyper- 
sonic limit. However, Kritsuk et al. (2007a) obtained Z2 ~ 0.76, 
whereas Z2 ~ 0.83 for our simulation, although the RMS Mach 
number is more than twice as high in the former case. Thus, 
our findings point to non-universal scaling properties of super- 
sonic turbulence. In addition to the RMS Mach number as major 
parameter, the intermittency of the compressible turbulence cas- 
cade might depend on the relative importance of compressive 
and solenoidal modes for the energy transfer among different 
scales. Indeed, the higher small-scale compressive ratio res ~ 0.5 
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in our simulation compared to about 0.3 in the simulation of 
Kritsuk et al. indicates a larger fraction of dilatational (rotation- 
free) velocity fluctuations although the RMS Mach number is 
lower Of course, this is a direct consequence of the applied 
forcing as far as large-scale velocity fluctuations are concerned. 
However, the divergence and the vorticity are mostly structural 
probes of small-scale velocity fluctuations (Frisch 1995). For 
this reason, the result for res suggests an influence of the forc- 
ing on the small-scale structure of turbulence. Remarkably, even 
in the incompressible limit, non-universal scaling at small scales 
has b een found in numerical simulations of wall-bounded shear 
flow jCasciola et al. 2007) and for turb ulence between counter- 
rotating cyhnders lAlexakis et alj|2005h . Particularly interesting 
is the spectral analysis presented in the latter work. It is shown 
that the turbulence energy flux at small scales is strongly influ- 
enced by triadic interactions with large scale modes and, thus, 
even distant scales cannot be separated. 

The one-parameter model ( |29] l features the codimension C 
as a measure of intermittency depending on the RMS Mach 
number, but the scaling parameter A - 2/3 stems from the 
incompressible turbulence cascade. However, it is not at all 
clear whether this scaling applies to supersonic turbulence. 
Consequently, we fitted the general two-parameter model (IZTT i 
subject to the constraint A < 1 to our data up to fourth order. 
The resulting parameters are A ^ 1.00 and p x 0.57 < 2/3. 
We also investigated the influence of including Z5 in the data set 
for the fit and found that the results did not change significantly 
(the value of C varied by less than 0.02). However, releasing 
the constraint A < 1, a value slightly greater than unity is ob- 
tained for A. Extending the fit to the scaling exponents up to 
sixth Oder, yields A » 1.05. Therefore, higher-order appears to 
be essential to constrain the A-parameter. Remarkably, A = 1 
follows from the scaling argume nt for the m ost singular dissipa- 
tive structures that was used by IShe & Levequd d 19941) . if it is 
assumed that the uniform dissipation time scale is proportional 
to £ as in Burgers turbulenceQ Therefore, the two-parameter fit 
implies shocks as the most singular dissipative structures. The 
/3-parameter indicates a higher degree of intermittency in com- 
parison to the Kolmogorov-Burgers model. The corresponding 
codimension C ^ 1.23 is in accordance with the notion that the 
codimension varies form C = 2 in the subsonic and C = 1 in the 
hypersonic limit. 

Whereas both the one-parameter model with C ~ 0.7 1 and 
the two-parameter model with C ~ 1 .23 agree excellently with 
the exponents Zi and Z2, the latter matches the higher-order 
statistics much better. However, the scaling exponents for p > 3 
are subject to greater uncertainty (see Figure [TSll. For this rea- 
son, we cannot C = 0.71. Interestingly, the sa me dichotomy 
with reg ard to A = 2/3 and A = 1 was found by iKritsuk et al.l 
(l2007bl) . albeit for different values of /3 and C. In this respect, it 
is interesting to investigate the dependence of the codimension 
C on Z2, which is related to the scaling index y via Z2 - 27/^3. 
Substituting y6 = 1 - A/C and inverting equation (|27T i for p -2, 
the functions C(2y/(j,) for A = 2/3 and A = 1, respectively, are 
obtained. As one can see form the plots in Figure[T6] the solution 
for A = 1 converges towards C = 1 as Z2 — > 1 . Thus, the Burgers 
scaling follows from this model for 7 = 0.5 and ^3 = 1. This 
corresponds to the interpretation that dissipation is completely 
dominated by shocks in the hypersonic limit. On the other hand, 
for A = 2/3, which includes the Kolmogorov-Burger model with 




Fig. 15. Fits of the relative scaling exponents Zp - ^p/^3 to sev- 
eral models. The error bars indicate the scatter of the numeri- 
cally calculated exponents, the dots specify the averages. The 
thin dashed lines show the scalings predicted by She & Levequa 
( 1994) and by the Kolmogorov-Burger models ( |28] ) of lBoldvrevI 
(2002), r especti vely. The general log-Poisson model ( |27] ) of 
Dubrulle d 19941) yields the two fits shown as thick solid lines. 
The fit function that closely matches the higher-order exponents 
corresponds to A - 1, whereas the other function was obtained 
with the constraint A = 2/3. 
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^ The dissipation time scale is given by the kinetic energy ~ u'^(£) 
divided by the scale-invariant mean dissipation. For Burgers turbulence, 



Fig. 16. Codimension C - A/(j6 - l)as function of the second- 
order scaling index Z2 = 2y/^3 for the log-Poisson models with 
A - 2/3 (lower line ) and A = 1 (upper line), respectively. The 
two dots correspond to the fits to our numerical results shown in 
Figure [15] The Kolmogorov-Burgers model with C = 1 is also 
marked by a dot. 



C(0.74) - 2, no real solution exists if Z2 becomes greater than 
about 0.89 (corresponding to 7 = 0.45). The minimal codimen- 
sion is slightly less than 0.7. Judging from the limiting behaviour 
of both models, it appears that the log-Poisson model with A = 1 
is more sensible for supersonic turbulence. 

Recently, measurements of higher-order scalings for diffuse 
molecular clouds in Po laris and Taurus have become available 
jHilv-Blant et alj|2008l) . The relative scaling exponents obtained 
from these measurements, Z2 ~ 0.69 . . .0.71, Z4 x 1.27 .. . 1.30 
and Z5 ^ 1.53 .. . 1.60, fall in between the incompressible She- 
Leveque model and the Kolmogorov-Burgers model. The value 
Z2 ~ 0.7 can only be accommodated within the A = 2/3 branch. 
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Fig. 17. Temporally averaged probability density function of the 
logarithmic mass density fluctuations, ln(p/po), and log-normal 
fit functions (thin dashed lines). 



for which the codimension C is close to 2 (Figure [T6b. This cor- 
responds to weakly compressible, shear-dominated turbulence, 
plainly in opposition to th e scen ario discussed in this article. 
However, iHilv-Blant et al I (l20()8b consider only three different 
clouds with scaling expon ents that are at the re very lower end 
of those reported by (Bru nt & Heyed l2002h . Consequently, a 
much broader sample of higher-order velocity statistics would 
be required in order to exclude the occurrence of compression- 
dominated turbulence in the ISM. 



4.2. Density Statistics 

Previous theoretical and numerical and studies favour log- 
normal statistics fo r the ma ss density of su personic turbu lence 
in isothermal gas ( Vazquez-Semadeni 199?; Padoan et al.'1997'; 
iPassot & Vazquez-Semadeni 1998; Padoan & Nordlund 2002) . 
This means that the distribution of ^ = ln(p/po) is Gaussian. 
Consistency with (p/po) - 1 requires the variance of 6 being re- 
lated to the mean via cr^ = -2(5) . However, as one can see in 
Figure [17] a normal distribution subject to this constraint does 
not properly fit the time-averaged 6-pdf from our simulation. 
The standard deviation of t he closest possible fit is os ~ 1.76. 
iPadoan & Nordlundl (|2002') proposed that erg can be parame- 
terised in terms of the RMS Mach number: 



(r2=ln(l+0.25ML)- 



(30) 



We have Mrms ~ 2.3 in the interval of time over which the pdf of 
6 is averaged. For this Mach number, the above formula yields 
cTj * 0.92. The corresponding pdf is also plotted for comparison 
in Figure [T7] Clearly, compressively driven turbulence produces 
a substantially broader range of density fluctuations th an pre- 
dicted by Padoan and Nordlund. iKritsuk et all (l2007bl) . on the 
other hand, found a narrower range in relation to the RMS Mach 
number 

Since we have sampled the (5-pdfs over many integral time 
scales, the shape of our time-averaged pdf is very likely gen- 
uine. In test simulations, we have also excluded the possibilities 
that it is purely a resolution effect or a consequence of the nearly 
isothermal approximation. A skewed distribution of density fluc- 
tuations has important consequences. In Figure[T8] the total mass 



of gas with density higher than a given threshold density, 

M(p) = (2L)M p'pdf(p')dp' =po(2L)M exp(5) pdf(5) d5', 
Jp Js 

(31) 

is plotted for the time-averaged probability density function 
pdf((J) from our simulation, closest log-normal pdf and the log- 
normal pdf with cr,5 X 0.92, respectively. It is palpable that cal- 
culating M{p) on the basis of the log-normal fit with cr^ ~ 1 .76 
in place of the numerical pdf implies enormously wrong mass 
fractions in the high-density peaks. On other hand, the predic- 
tion of M(p) based on the Padoan-Nordlund relation ( |30] | under- 
estimates the total mass of highly compressed gas by about one 
order of magnitude. For p > 100, a numerical cutoff can be dis- 
cerned, as the mass function corresponding to the numerical pdf 
plunges towards zero in the logarithmically scaled plot (b). The 
reason is that density fluctuations cannot become arbitrarily high 
due to the spatial discretisation and, because of the discretisation 
in time, the most intermittent events are too rare to show up in 
the statistics. The mass range affected by the numerical cutoff is 
M(p) < 10"^. We found that this range largely overlaps with the 
mass range for the computation of clum p mass spectra following 
the prescription bv lPadoan et al.l ( l2007h . The high sensitivity of 
mass spectra on numerical resolution has already been noted by 
Hennebelle & Audit (2007). Consequently, it is mandatory to go 
to substantially higher resolutions which, in turn, necessitates 
the application of adaptive mesh refinement. Unlike the AMR 
simulation at hand, extreme density peaks have to be followed 
to very high resolution and, at the same time, turbulent flow has 
to be treated self-consistently in less resolved regions of lower 
gas density. 

A justification for log-normal pdfs of the mass 
density based on analyti cal argu ments was given by 
'Passot & Vazquez-SemadeM (Il998h for one-dimensional 
isothermal gas dynamics without external forces. These argu- 
ments also apply to fully developed turbulence in the strict 
sense, i. e., three-dimensional turbulence that evolves without 
imposed constraints such as boundaries, external forces or 
viscosity. To what extent can this be carried over to turbulence 
in a periodic box with stochastic forcing? Setting the pressure 
P - CqP, where cq is the isothermal speed of sound, the pressure 
gradient divided by the mass density is given by 



1 
P 



VP = civs. 



(32) 



Thus, the balance laws ([T]i and (|2]i can be written in the form 



—6 = -d, 
Dt 

D 

Dt 



u = -clV6 + f. 



(33) 
(34) 



The first equation means that infinitesimal Lagrangian changes 
of the density are given by -d dt. It is argued that density fluc- 
tuations are built up in a hierarchical process, meaning that 6 
evolves as a random process, where infinitesimal increments 
6+dt (density enhancements) and decrements 6-dt (density re- 
ductions) are added with equal probability independent of the 
value of 6. The central limit theorem then implies a Gaussian 
distribution of 6. This is the just the Wiener process that was 
mentioned in Section|2] These arguments also apply to the three- 
dimensional case. However, the action of an external force can 
cause deviations from the log-normal statistics. Let us consider 
equation ( fTTT i for the rate of compression in the isothermal case. 
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Fig. 18. Integrated mass M{p) in gas of density higher than p in linear scaling (a) and logarithmic scaling (b). The thick solid line 
results from the time-averaged pdf obtained from our simulation. The closest log-normal distribution and the log-normal distribution 
implied by the Padoan-Nordlund relation [30] yield the dashed cures. While the former greatly overestimates the mass in the density 
peaks, the latter underestimates the fraction of over-dense gas. 



In regions, in which V ■ / < 0, the force tends to increase the rate 
of compression. The equation for the density fluctuations ( l33T l 
implies that the corresponding fractional change of -ddt re- 
sults in a stronger increment 6+dt or a smaller decrement 5-dt, 
i. e., a force with positive convergence locally supports com- 
pression and weakens rarefaction, as one would expect. In the 
case V ■ / > 0, the force has the opposite effect. If the force 
field is stochastic, both effects will occur with equal probabil- 
ity at any time at any position (by the very construction of the 
force field). Hence, the net effect depends on the periods of time 
a particular fluid parcel is contracting {-d > 0) or expanding 
{-d < 0). Apart from that, the effect will weaken toward smaller 
scales and higher densities (because of scaling considerations). 
For a detailed analysis, Lagrangian density statistics would be 
required. From the existing simulation, we have only Eulerian 
statistics. Nevertheless, our considerations show that, in general, 
forces with a strong dilatational component may cause devia- 
tions form a log-normal distribution of the density fluctuations 
in isothermal gas. This is what we observe in our simulation. 

The spectral properties of 6 in our simulation indicate an in- 
termittent self-similar structure, where the scaling appears to be 
closely linked to the scaling of the velocity fluctuations. The 
spectral index jSs ~ 1.84 implies an exponent ys ~ 0.42 of 
the RMS logarithmic density fluctuation. The corresponding in- 
dex of the RMS velocity fluctuations is y a; 0.43. Similar scal- 
ing properties of the velocity and density fluctuations are im- 
plied by equation ( |33] ) provided that / only acts on large scales. 
Therefore, self-similarity (related to scaling) appears to be a 
more robust property of turbulent density fluctuations than log- 
normality (a probabilistic concept), contrary to the explanation 
of log-normal statistics on grounds of self-similarity in the liter- 
ature. 



5. Conclusion 

We performed simulations of supersonic turbulence in isother- 
mal gas driven by a mostly compressive stochastic force field. 
The analysis of statistical properties of the velocity and mass 
density revealed several remarkable results: 



1 . The turbulence energy spectrum function is close to the spec- 
trum of Burgers turbulence in the initial phase of the flow 
evolution and subsequently approaches a flatter spectrum 
with an index [3 k 1.9. Within the error bars, a compara- 
ble index j3s ~ 1.8 is found for the spectrum function of 6^, 
where 6 - ln(p/po) is the logarithmic density fluctuation. 

2. In the statistically stationary state, the scaling exponent of 
the RMS velocity fluctuation y = ^2/2 ~ 0.43, which is com- 
parable to typical exponents inferred from measurements of 
the velocity dispersion in molecular clouds. 

3. The relative scaling exponents Zp - up to order p - 5 
deviate from the Kolmogorov-Burgers model for supersonic 
turbulence. Varying only the codimension C of dissipative 
structures, the best fit to our data implies C ~ 0.71, whereas 
C = 1 in the Kolmogorov-Burgers model. Fitting the general 
log-Poisson hierarchical structure model, we find C ~ 1.23 
and the scaling index A x 1.0. This can be interpreted as 
shocks being the most singular dissipative structures. In the 
Kolmogorov-Burgers model, A - 2/3 is assumed. 

4. Compared to other simulations, our analysis point toward 
non-universal scaling properties, i. e., the velocity statistics 
cannot be accommodated within a one-parameter family of 
models. There appears to be an additional degree of freedom 
related to the large-scale forcing. 

5. The time-averaged probability density function of 6 is not 
consistent with a log-normal distribution of the mass density. 
We find a skewed distribution t hat is broader than the log- 
normal distribution predicted bv'P adoan & Nordlundl (|2002|) 
on the basis of the RMS Mach number Large deviations of 
the mass fraction in over-dense gas compared to log-normal 
distributions are implied. 

The assumption of a log-normal mass density distribu- 
tion has been used for analytical predictions of the rate of 
turbulence-regulated star formation jP adoan & Nordlund 2002; 
iKmmholz & McKeell2005l; iHennebeUe & Chabrien i2008). With 
the distribution of the mass density found in our simulation, 
the resulting star formation rate would differ substantially form 
the log-normal case because the mass converted into stars by 
gravitational collapse mainly depends on the high-density tail 
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of the probability distribution function. However, this result is 
not conclusive yet, because a wider range of parameters has to 
be studied, self-gravity and magnetic fields have to be included, 
and, most importantly, the dependence of the mass density dis- 
tribution and clump mass spectra on the numerical resolution 
has to be clarified. In part, this is the subject of ongoing work, 
where turbulence simulations with purely solenoidal and com- 
pressive forcing are systematically compared for numerical res - 
olutions up to 1024^ dSchmidt et al.l2008HFederrath et aljioosh . 
Moreover, this complementary study is essential to corroborate 
the non-universal velocity scaling properties discussed in this ar- 
ticle. 

The observational results of iBrunt & Heveil (|2002|) pomt to- 
wards a scaling exponent a that is significantly greater than 
0.5 in many instances. This is hard to explain on the ba- 
sis of supersonic turbulence in isothermal gas. StifFer scal- 
ing exp onents might be caused, for instance , by self-gravity 
(Biglari & Diamondlll989h . On the other hand. iHilv-Blant et all 
(j2008) found scaling properties that are indicative of nearly in- 
compressible, shear-dominated turbulence in some molecular 
clouds. Possibly, molecular clouds are a rather inhomogeneous 
class of objects. From that point of view, the situation is rem- 
iniscent of the observation of type la supernovae, which were 
recognised to be much more diverse than it had appeared in 
the beginning. To understand how the mass density distribution 
and the two-point statistics of the turbulent velocity at various 
scales in the interstellar medium comes about, it is vital to in- 
clude additional physics such as self-gravity, radiative cooling 
and magnetic fields also in numerical simulations of turbulence 
with stochastic forcing. To this e nd, the appl ication of adaptive 
methods will be indispensable. lOff'ner et al.l (j2008) succeeded 
in performing simulations of self-gravitating isothermal turbu- 
lence with AMR. Radiative post-processing of their numerical 
data revealed discrepancies with observed core line widths both 
for driven and for decaying turbulence. This suggests that not 
only gravity but also thermal processes are important for the 
dynamics of molecular cloud turbulence. With the AMR sim- 
ulation presented in this article, we have attempted a first step 
toward a methodology tha t is based on the dynamics of turbu- 
lent flows. Ilapichino et"aLl (12008.) have already presented a multi- 
level AMR application using these techniques. The development 
is far from being complete yet, but it appears to be promising for 
a self-consistent treatment of turbulence in the complex scenar- 
ios encountered in star-forming clouds. 
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